Here I will explore how DEPOSITED AND SUSPENDED sediment affects the MORTALITY of tropical, scleractinean coral ADULTS. The database to be used is comprised of data extracted from studies deemed relevant during the systematic literature review process.

The Dataset

Deposited Sediment

Specifically, we will explore the results of 827 datapoints from 28 studies conducted in 3 oceans on 87 species within 46 genera.

##     n
## 1 827
##      Ref   n
## 1  DS08a 162
## 2  DS08b 113
## 3   DS71 112
## 4   DS04  70
## 5   DS68  48
## 6   DS46  44
## 7   DS29  40
## 8   DS30  36
## 9   DS17  24
## 10  DS45  24
## 11  DS16  19
## 12  DS18  16
## 13  DS27  16
## 14  DS32  16
## 15  DS07  13
## 16  DS03  12
## 17  DS43  12
## 18  DS05  10
## 19  DS33   8
## 20  DS12   7
## 21  DS26   6
## 22  DS11   4
## 23  DS21   4
## 24  DS28   4
## 25  DS25   2
## 26  DS36   2
## 27  DS69   2
## 28  DS23   1
##                                       Ref_name   n
## 1                Hodgson (1989) PhD Chapter IV 275
## 2  HDR EOC and CSA Ocean Services, Inc. (2014) 112
## 3                      Duckworth et al. (2017)  70
## 4                         Flores et al. (2012)  48
## 5                        Stafford-Smith (1993)  44
## 6          Stafford-Smith (1990) PhD Chapter 4  40
## 7                       Stafford-Smith (1992)   36
## 8                                Piniak (2007)  24
## 9                                Rogers (1983)  24
## 10                Philipp and Fabricius (2003)  19
## 11                     Piniak and Brown (2008)  16
## 12                Sofonia (2006) PhD Chapter 4  16
## 13                       Stewart et al. (2006)  16
## 14                             Hodgson (1990b)  13
## 15              Bessell-Browne et al. (2017a)   12
## 16                   Peters and Pilson (1985)   12
## 17                                Hodel (2007)  10
## 18                  Vargas-Angel et al. (2006)   8
## 19                        Loiola et al. (2013)   7
## 20                Sofonia (2006) PhD Chapter 3   6
## 21                        Lirman et al. (2008)   4
## 22                              Rogers (1979)    4
## 23                  Sofonia and Anthony (2008)   4
## 24                           Gil et al. (2016)   2
## 25                  Shore-Maggio et al. (2018)   2
## 26                          Zill et al. (2017)   2
## 27                                Selim (2007)   1
##      Ocean   n
## 1  Pacific 757
## 2 Atlantic  69
## 3   Indian   1
## Selecting by n
##                           Gsp  n
## 1          Acropora_millepora 38
## 2            Leptoria_phrygia 42
## 3  Montipora_aequituberculata 26
## 4              Oxypora_glabra 23
## 5               Pavona_cactus 29
## 6          Porites_cylindrica 30
## 7        Porites_lobata/lutea 56
## 8               Porites_lutea 28
## 9                 Porites_rus 30
## 10      Turbinaria_reniformis 24
## Selecting by n
##    Updated_Genus   n
## 1        Porites 152
## 2      Montipora 137
## 3       Acropora  88
## 4       Leptoria  42
## 5         Pavona  37
## 6     Turbinaria  33
## 7        Oxypora  31
## 8    Pocillopora  29
## 9     Echinopora  23
## 10      Mycedium  21

Suspended Sediment

Specifically, we will explore the results of 272 datapoints from 6 studies conducted in 3 oceans on 19 species within 14 genera (shown below only the top 10 species and genera, by number of datapoints).

##     n
## 1 272
##     Ref   n
## 1  SS06 150
## 2  SS08  48
## 3 SS22c  28
## 4  SS27  16
## 5 SS22a  10
## 6  SS04   9
## 7 SS22b   8
## 8  SS25   3
##                     Ref_name   n
## 1         Browne et al. 2015 150
## 2         Flores et al. 2012  48
## 3                  Rice 1984  46
## 4        Anthony et al. 2007  16
## 5 Bessell-Browne et al. 2017   9
## 6          Te 2001 Chapter 6   3
##      Ocean   n
## 1   Indian 150
## 2  Pacific  76
## 3 Atlantic  46
## Selecting by n
##                           Gsp  n
## 1           Merulina_ampliata 50
## 2         Pachyseris_speciosa 50
## 3          Platygyra_sinensis 50
## 4  Montipora_aequituberculata 30
## 5          Acropora_millepora 21
## 6         Acropora_intermedia 16
## 7         Cladocora_arbuscula  8
## 8           Manicina_areolata  8
## 9          Solenastrea_hyades  8
## 10        Siderastrea_radians  6
## Selecting by n
##     Updated_Genus  n
## 1        Merulina 50
## 2      Pachyseris 50
## 3       Platygyra 50
## 4        Acropora 37
## 5       Montipora 36
## 6       Cladocora  8
## 7        Manicina  8
## 8     Solenastrea  8
## 9     Siderastrea  6
## 10     Isophyllia  4
## 11     Phyllangia  4
## 12       Scolymia  4
## 13 Stephanocoenia  4

Exploratory plots

First let’s explore all data from all species for which there is binary data about the presence of ‘small necroses’, ‘large necroses’, and ‘total mortality’ as a result of exposure to deposited sediment.

DEFINITIONS:
  • ‘Small necroses’ are defined as areas of tissue mortality that encompass between 0% and 50% of the coral’s surface area.
  • ‘Large necroses’ are defined as areas of tissue mortality that encompass greater than or equal to 50% of the coral’s surface area, but less than 100% (total mortality).
  • ‘Total mortality’ is defined as 100% of coral’s surface area experiencing necrosis. In other words, it is the death of an entire coral entity (e.g., colony, branch, fragment, other experimental unit).

Threshold Analyses with Binary Data

Now let’s calculate the thresholds for each category of mortality, based on the binary data explored above.

DEFINITIONS:
  • ‘LOAEL’, or the ‘lowest observed adverse effect level’ is defined as the lowest dose at which there was an observed adverse effect.
  • ‘NOAEL’, or the ‘no observed adverse effect level’ is defined as the highest dose at which there was NOT an observed adverse effect.
  • ‘Adverse effect’ is defined here as any response of a coral individual, colony, or experimental treatment group that may negatively affect the coral’s fitness and/or survival. These adverse effects may include sublethal physiological changes (e.g., significantly reduced growth or photosynthetic rates when compared to coral in ambient/control conditions), bleaching/tissue loss, and mortality. This definition of adverse effect is independent of its magnitude, so while the affect may have potentially reduced fitness, its magnitude may be sufficiently small that the fitness effect is not measureable.

DEPOSITED SEDIMENT

Small Necroses

NOAEL/LOAEL defined by exposure concentration

##   LOAEL
## 1   4.9
##   NOAEL
## 1   4.4

NOAEL/LOAEL defined by exposure duration

##   LOAEL
## 1 0.917
##   NOAEL
## 1 0.917

Large Necroses

NOAEL/LOAEL defined by exposure concentration

##   LOAEL
## 1  20.8
##   NOAEL
## 1  20.8

NOAEL/LOAEL defined by exposure duration

##   LOAEL
## 1     3
##   NOAEL
## 1     3

Total mortality

NOAEL/LOAEL defined by exposure concentration

##   LOAEL
## 1  20.8
##   NOAEL
## 1  20.8

NOAEL/LOAEL defined by exposure duration

##   LOAEL
## 1     1
##   NOAEL
## 1     1

SUSPENDED SEDIMENT

Small Necroses

NOAEL/LOAEL defined by exposure concentration

##   LOAEL
## 1  3.22
##   NOAEL
## 1  3.22

NOAEL/LOAEL defined by exposure duration

##   LOAEL
## 1    14
##   NOAEL
## 1    14

Large Necroses

NOAEL/LOAEL defined by exposure concentration

##   LOAEL
## 1  29.1
##   NOAEL
## 1  29.1

NOAEL/LOAEL defined by exposure duration

##   LOAEL
## 1    84
##   NOAEL
## 1    84

Total mortality

NOAEL/LOAEL defined by exposure concentration

##   LOAEL
## 1  29.1
##   NOAEL
## 1  29.1

NOAEL/LOAEL defined by exposure duration

##   LOAEL
## 1    84
##   NOAEL
## 1    84

Logistic Regression Model Fitting

DEPOSITED SEDIMENT

Small Necroses

##     n
## 1 602
## Data: smallDS2
## Models:
## modDS6: Binary_small_necroses ~ Sed_level_stand_mg * Sed_exposure_days + 
## modDS6:     (1 | Gsp)
## modDS9: Binary_small_necroses ~ Sed_level_stand_mg * Sed_exposure_days + 
## modDS9:     (1 | Ref)
## modDS12: Binary_small_necroses ~ Sed_level_stand_mg * Sed_exposure_days + 
## modDS12:     (1 | Gsp) + (1 | Ref)
##         npar    AIC    BIC  logLik deviance  Chisq Df Pr(>Chisq)    
## modDS6     5 660.76 682.77 -325.38   650.76                         
## modDS9     5 702.81 724.81 -346.40   692.81  0.000  0          1    
## modDS12    6 652.31 678.71 -320.15   640.31 52.496  1   4.31e-13 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## Data: smallDS2
## Models:
## modDS10: Binary_small_necroses ~ Sed_level_stand_mg + (1 | Gsp) + (1 | 
## modDS10:     Ref)
## modDS11: Binary_small_necroses ~ Sed_level_stand_mg + Sed_exposure_days + 
## modDS11:     (1 | Gsp) + (1 | Ref)
## modDS12: Binary_small_necroses ~ Sed_level_stand_mg * Sed_exposure_days + 
## modDS12:     (1 | Gsp) + (1 | Ref)
##         npar    AIC    BIC  logLik deviance  Chisq Df Pr(>Chisq)   
## modDS10    4 657.21 674.81 -324.61   649.21                        
## modDS11    5 650.97 672.97 -320.48   640.97 8.2432  1    0.00409 **
## modDS12    6 652.31 678.71 -320.15   640.31 0.6604  1    0.41641   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: binomial  ( logit )
## Formula: Binary_small_necroses ~ Sed_level_stand_mg + Sed_exposure_days +  
##     (1 | Gsp) + (1 | Ref)
##    Data: smallDS2
## 
##      AIC      BIC   logLik deviance df.resid 
##    651.0    673.0   -320.5    641.0      597 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.1315 -0.6021 -0.1801  0.5730  3.0137 
## 
## Random effects:
##  Groups Name        Variance Std.Dev.
##  Gsp    (Intercept) 2.871    1.694   
##  Ref    (Intercept) 2.404    1.551   
## Number of obs: 602, groups:  Gsp, 75; Ref, 21
## 
## Fixed effects:
##                     Estimate Std. Error z value Pr(>|z|)    
## (Intercept)        -1.975421   0.562799  -3.510 0.000448 ***
## Sed_level_stand_mg  0.002889   0.001178   2.453 0.014167 *  
## Sed_exposure_days   0.013832   0.004908   2.818 0.004828 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## convergence code: 0
## Model is nearly unidentifiable: very large eigenvalue
##  - Rescale variables?
## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: binomial  ( logit )
## Formula: Binary_small_necroses ~ log10_conc + log10_dur + (1 | Gsp) +  
##     (1 | Ref)
##    Data: smallDS2
## 
##      AIC      BIC   logLik deviance df.resid 
##    639.8    661.8   -314.9    629.8      597 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.3280 -0.5459 -0.1809  0.5600  3.1754 
## 
## Random effects:
##  Groups Name        Variance Std.Dev.
##  Gsp    (Intercept) 3.219    1.794   
##  Ref    (Intercept) 3.934    1.983   
## Number of obs: 602, groups:  Gsp, 75; Ref, 21
## 
## Fixed effects:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  -5.4867     1.1307  -4.852 1.22e-06 ***
## log10_conc    1.3053     0.3869   3.374 0.000742 ***
## log10_dur     1.5651     0.4175   3.749 0.000178 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Effect estimates

## [1] 0.8272425
##    
## p     0   1
##   0 344  61
##   1  43 154
## Setting levels: control = 0, case = 1
## Setting direction: controls < cases
## Area under the curve: 0.9049
##                    R2m       R2c
## theoretical 0.07964959 0.7100653
## delta       0.07245285 0.6459074

## $Gsp

## 
## $Ref

##  Groups Name        Std.Dev.
##  Gsp    (Intercept) 1.7942  
##  Ref    (Intercept) 1.9834
## [1] 7.267637

##                     Est           LL          UL
## (Intercept) 0.004141512 0.0004515197  0.03798754
## log10_conc  3.688788053 1.7279975250  7.87452360
## log10_dur   4.783110717 2.1103679520 10.84083376

For every 10-fold increase in exposure concentration of deposited sediment, the odds of developing small necroses (<50% of coral tissue surface area) increase by 3.7 times (95% CI 1.7, 7.9; GLMM z = 3.374, p = 0.0007), after accounting for exposure duration and variability among studies and species. For every 10-fold increase in exposure duration of deposited sediment, the odds of developing small necroses increase by 4.8 times (95% CI 2.1, 10.8; GLMM z = 3.749, p = 0.0002), after accounting for exposure concentration and variability among studies and species.

Prediction figure

#Overlaying glmm results on figure, but prediction line is jagged...
pred_modDS11_log <- predict(modDS11_log, newdata=smallDS2, type="response")
smallDS3 <- cbind(smallDS2, pred_modDS11_log)

smallDS_plot <- ggplot(data = smallDS3) +
    geom_point(mapping = aes(
      x = Sed_level_stand_mg, 
      y = Binary_small_necroses)) +
    labs(x = expression("Sediment exposure concentration (mg/cm"^"2"*"/day)"), 
         y = "Small necroses due to sediment exposure?") +
    scale_x_log10(limits=c(1,1000),breaks=c(0.01,0.1,1,10,100,1000),
                  label=c("0.01","0.1","1","10","100","1000")) +
    geom_line(aes(x = Sed_level_stand_mg, y = pred_modDS11_log), inherit.aes=FALSE) +
    theme_bw()

smallDS_plot

That plot is confusing to interpret. Let’s see if I can plot the average marginal probability, i.e., the average change in probability of the outcome across the range of the predictor of interest. This is described in some detail at the following, useful website: https://stats.idre.ucla.edu/r/dae/mixed-effects-logistic-regression/

By sediment exposure concentration
jvalues <- with(smallDS2, seq(from = min(log10_conc), to = max(log10_conc), length.out = 100))

# calculate predicted probabilities and store in a list
pp <- lapply(jvalues, function(j) {
    smallDS2$log10_conc <- j
    predict(modDS11_log, newdata = smallDS2, type = "response")
})

# average marginal predicted probability across a few different sediment levels
#sapply(pp[c(1, 20, 40, 60, 80, 100)], mean)
## 
# get the means with lower and upper quartiles
plotdat <- t(sapply(pp, function(x) {
    c(M = mean(x), med = median(x), quantile(x, c(0.25, 0.75)))
}))

# add in log10_conc values and convert to data frame
plotdat <- as.data.frame(cbind(plotdat, jvalues))
plotdat <- plotdat %>% mutate(Sed_level_linear=10^jvalues)

# better names and show the first few rows
colnames(plotdat) <- c("PredictedProbability", "MedianProbability", 
                       "LowerQuantile", "UpperQuantile", "log10_conc", "Sed_conc")
#head(plotdat)

# plot average marginal predicted probabilities
ggplot(plotdat, aes(x = log10_conc)) + 
  geom_line(aes(y = PredictedProbability), size = 2) +
  geom_line(aes(y = MedianProbability), size = 0.5) +
  ylim(c(0, 1))

ggplot(plotdat, aes(x = log10_conc)) + 
  geom_ribbon(aes(ymin = LowerQuantile, ymax = UpperQuantile), alpha = 0.15) +
  geom_line(aes(y = PredictedProbability), size = 2) +
  geom_line(aes(y = MedianProbability), size = 0.5) + 
  ylim(c(0, 1))

ggplot(plotdat, aes(x = Sed_conc, y = PredictedProbability)) + 
  geom_ribbon(aes(ymin = LowerQuantile, ymax = UpperQuantile), alpha = 0.15) +
  geom_line(aes(y = PredictedProbability), size = 2) +
  geom_line(aes(y = MedianProbability), size = 0.5) + 
  ylim(c(0, 1)) +
  scale_x_log10()

#Overlaying average marginal predicted probabilities on figure
#by Ref
smallDS_plot2 <- ggplot() +
    geom_point(data = smallDS2, mapping = aes(
      x = Sed_level_stand_mg, 
      y = Binary_small_necroses,
      color = Ref)) +
    labs(x = expression("Sediment exposure concentration (mg/cm"^"2"*"/day)"), 
         y = "Predicted probability of small\nnecroses due to sediment exposure",
         color = "\nBinary Data by Study",
         linetype = "Predicted\nProbability") +
    scale_x_log10(limits=c(1,1000), breaks=c(0.01,0.1,1,10,100,1000,loaelDS1),
                  label=c("0.01","0.1","1","10","100","1000",round(loaelDS1,digits=1))) +
    geom_ribbon(data = plotdat, aes(x = Sed_conc, y = PredictedProbability, 
                                    ymin = LowerQuantile, ymax = UpperQuantile), 
                alpha = 0.15) +
    geom_line(data = plotdat, aes(x = Sed_conc, y = PredictedProbability, 
                                  linetype = "twodash")) +
    geom_line(data = plotdat, aes(x = Sed_conc, y = MedianProbability, 
                                  linetype = "solid")) +
    geom_vline(xintercept=loaelDS1, linetype="dashed", color = "red") +
    theme_bw() +
    scale_linetype_manual(values=c("twodash", "solid"), labels = c("Median","Mean"))

smallDS_plot2

By sediment exposure duration
jvalues <- with(smallDS2, seq(from = min(log10_dur), to = max(log10_dur), length.out = 100))

# calculate predicted probabilities and store in a list
pp <- lapply(jvalues, function(j) {
    smallDS2$log10_dur <- j
    predict(modDS11_log, newdata = smallDS2, type = "response")
})

# average marginal predicted probability across a few different sediment levels
#sapply(pp[c(1, 20, 40, 60, 80, 100)], mean)
## [1] 0.1797186 0.1960244 0.2142042 0.2333993 0.2535694 0.2746604
# get the means with lower and upper quartiles
plotdat <- t(sapply(pp, function(x) {
    c(M = mean(x), med = median(x), quantile(x, c(0.25, 0.75)))
}))

# add in Sed_exposure_days values and convert to data frame
plotdat <- as.data.frame(cbind(plotdat, jvalues))
plotdat <- plotdat %>% mutate(Sed_dur_linear=10^jvalues)

# better names and show the first few rows
colnames(plotdat) <- c("PredictedProbability", "MedianProbability", 
                       "LowerQuantile", "UpperQuantile", "log10_dur","Sed_dur")
#head(plotdat)

# plot average marginal predicted probabilities
ggplot(plotdat, aes(x = log10_dur)) + 
  geom_line(aes(y = PredictedProbability), size = 2) +
  geom_line(aes(y = MedianProbability), size = 0.5) +
  ylim(c(0, 1))

ggplot(plotdat, aes(x = log10_dur)) + 
  geom_ribbon(aes(ymin = LowerQuantile, ymax = UpperQuantile), alpha = 0.15) +
  geom_line(aes(y = PredictedProbability), size = 2) +
  geom_line(aes(y = MedianProbability), size = 0.5) + 
  ylim(c(0, 1))

ggplot(plotdat, aes(x = Sed_dur, y = PredictedProbability)) + 
  geom_ribbon(aes(ymin = LowerQuantile, ymax = UpperQuantile), alpha = 0.15) +
  geom_line(aes(y = PredictedProbability), size = 2) +
  geom_line(aes(y = MedianProbability), size = 0.5) + 
  ylim(c(0, 1)) +
  scale_x_log10()

#Overlaying average marginal predicted probabilities on figure
#by Ref
smallDS_plot3 <- ggplot() +
    geom_point(data = smallDS2, mapping = aes(
      x = Sed_exposure_days, 
      y = Binary_small_necroses,
      color = Ref)) +
    labs(x = "Sediment exposure duration (days)", 
         y = "Predicted probability of small\nnecroses to sediment exposure",
         color = "\nBinary Data by Study",
         linetype = "Predicted\nProbability") +
    scale_x_log10(limits=c(0.1,200), breaks=c(0.01,0.1,1,10,100,1000),
                  label=c("0.01","0.1","1","10","100","1000")) +
    geom_ribbon(data = plotdat, aes(x = Sed_dur, y = PredictedProbability, 
                                    ymin = LowerQuantile, ymax = UpperQuantile), 
                alpha = 0.15) +
    geom_line(data = plotdat, aes(x = Sed_dur, y = PredictedProbability, 
                                  linetype = "twodash")) +
    geom_line(data = plotdat, aes(x = Sed_dur, y = MedianProbability, 
                                  linetype = "solid")) +
    theme_bw() +
    scale_linetype_manual(values=c("twodash", "solid"), labels = c("Median","Mean"))

smallDS_plot3

Large Necroses

##     n
## 1 522
## boundary (singular) fit: see ?isSingular
## Data: largeDS2
## Models:
## modDSb6: Binary_large_necroses ~ Sed_level_stand_mg * Sed_exposure_days + 
## modDSb6:     (1 | Gsp)
## modDSb9: Binary_large_necroses ~ Sed_level_stand_mg * Sed_exposure_days + 
## modDSb9:     (1 | Ref)
## modDSb12: Binary_large_necroses ~ Sed_level_stand_mg * Sed_exposure_days + 
## modDSb12:     (1 | Gsp) + (1 | Ref)
##          npar    AIC    BIC  logLik deviance  Chisq Df Pr(>Chisq)    
## modDSb6     5 346.09 367.38 -168.04   336.09                         
## modDSb9     5 392.60 413.88 -191.30   382.60  0.000  0          1    
## modDSb12    6 345.25 370.80 -166.62   333.25 49.346  1  2.146e-12 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## Data: largeDS2
## Models:
## modDSb6: Binary_large_necroses ~ Sed_level_stand_mg * Sed_exposure_days + 
## modDSb6:     (1 | Gsp)
## modDSb12: Binary_large_necroses ~ Sed_level_stand_mg * Sed_exposure_days + 
## modDSb12:     (1 | Gsp) + (1 | Ref)
##          npar    AIC    BIC  logLik deviance  Chisq Df Pr(>Chisq)  
## modDSb6     5 346.09 367.38 -168.04   336.09                       
## modDSb12    6 345.25 370.80 -166.62   333.25 2.8366  1    0.09214 .
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## Data: largeDS2
## Models:
## modDSb4: Binary_large_necroses ~ Sed_level_stand_mg + (1 | Gsp)
## modDSb5: Binary_large_necroses ~ Sed_level_stand_mg + Sed_exposure_days + 
## modDSb5:     (1 | Gsp)
## modDSb6: Binary_large_necroses ~ Sed_level_stand_mg * Sed_exposure_days + 
## modDSb6:     (1 | Gsp)
##         npar    AIC    BIC  logLik deviance  Chisq Df Pr(>Chisq)   
## modDSb4    3 350.54 363.31 -172.27   344.54                        
## modDSb5    4 345.35 362.38 -168.68   337.35 7.1909  1   0.007327 **
## modDSb6    5 346.09 367.38 -168.04   336.09 1.2629  1   0.261106   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## Data: largeDS2
## Models:
## modDSb10: Binary_large_necroses ~ Sed_level_stand_mg + (1 | Gsp) + (1 | 
## modDSb10:     Ref)
## modDSb11: Binary_large_necroses ~ Sed_level_stand_mg + Sed_exposure_days + 
## modDSb11:     (1 | Gsp) + (1 | Ref)
## modDSb12: Binary_large_necroses ~ Sed_level_stand_mg * Sed_exposure_days + 
## modDSb12:     (1 | Gsp) + (1 | Ref)
##          npar    AIC    BIC  logLik deviance  Chisq Df Pr(>Chisq)   
## modDSb10    4 352.54 369.57 -172.27   344.54                        
## modDSb11    5 343.95 365.24 -166.98   333.95 10.589  1   0.001137 **
## modDSb12    6 345.25 370.80 -166.62   333.25  0.701  1   0.402432   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## Data: largeDS2
## Models:
## modDSb5: Binary_large_necroses ~ Sed_level_stand_mg + Sed_exposure_days + 
## modDSb5:     (1 | Gsp)
## modDSb11: Binary_large_necroses ~ Sed_level_stand_mg + Sed_exposure_days + 
## modDSb11:     (1 | Gsp) + (1 | Ref)
##          npar    AIC    BIC  logLik deviance  Chisq Df Pr(>Chisq)  
## modDSb5     4 345.35 362.38 -168.68   337.35                       
## modDSb11    5 343.95 365.24 -166.98   333.95 3.3985  1    0.06526 .
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: binomial  ( logit )
## Formula: Binary_large_necroses ~ Sed_level_stand_mg + Sed_exposure_days +  
##     (1 | Gsp) + (1 | Ref)
##    Data: largeDS2
## 
##      AIC      BIC   logLik deviance df.resid 
##    344.0    365.2   -167.0    334.0      517 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.9660 -0.3198 -0.0383 -0.0265  4.5830 
## 
## Random effects:
##  Groups Name        Variance Std.Dev.
##  Gsp    (Intercept) 28.3671  5.3261  
##  Ref    (Intercept)  0.9538  0.9766  
## Number of obs: 522, groups:  Gsp, 74; Ref, 18
## 
## Fixed effects:
##                     Estimate Std. Error z value Pr(>|z|)    
## (Intercept)        -7.702214   2.000071  -3.851 0.000118 ***
## Sed_level_stand_mg  0.003963   0.001377   2.878 0.003998 ** 
## Sed_exposure_days   0.026231   0.009663   2.715 0.006635 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## convergence code: 0
## Model is nearly unidentifiable: very large eigenvalue
##  - Rescale variables?
## Model is nearly unidentifiable: large eigenvalue ratio
##  - Rescale variables?
## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: binomial  ( logit )
## Formula: Binary_large_necroses ~ log10_conc + log10_dur + (1 | Gsp) +  
##     (1 | Ref)
##    Data: largeDS2
## 
##      AIC      BIC   logLik deviance df.resid 
##    322.9    344.2   -156.5    312.9      517 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.7377 -0.1869 -0.0426 -0.0201  4.1335 
## 
## Random effects:
##  Groups Name        Variance Std.Dev.
##  Gsp    (Intercept) 27.088   5.205   
##  Ref    (Intercept)  1.146   1.071   
## Number of obs: 522, groups:  Gsp, 74; Ref, 18
## 
## Fixed effects:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept) -14.0241     2.7141  -5.167 2.38e-07 ***
## log10_conc    2.5721     0.5767   4.460 8.20e-06 ***
## log10_dur     2.7644     0.8983   3.077  0.00209 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Effect estimates

## [1] 0.9099617
##    
## p     0   1
##   0 434  34
##   1  13  41
## Setting levels: control = 0, case = 1
## Setting direction: controls < cases
## Area under the curve: 0.9494
##                    R2m       R2c
## theoretical 0.09338887 0.9053851
## delta       0.09041040 0.8765094

## $Gsp

## 
## $Ref

##  Groups Name        Std.Dev.
##  Gsp    (Intercept) 5.2046  
##  Ref    (Intercept) 1.0707
## [1] 2.917504

##                      Est           LL           UL
## (Intercept) 8.117388e-07 3.973205e-09 1.658409e-04
## log10_conc  1.309276e+01 4.227962e+00 4.054444e+01
## log10_dur   1.586911e+01 2.728155e+00 9.230725e+01

For every 10-fold increase in exposure concentration of deposited sediment, the odds of any tissue mortality increase by 13.1 times (95% CI 4.2, 40.5; GLMM z = 4.460, p < 0.00001), while holding exposure duration constant and accounting for variability among studies and species. For every 10-fold change in exposure duration, the odds of any tissue mortality increase by 15.9 times (95% CI 2.7, 9.2; GLMM z = 3.077, p = 0.002), while holding exposure concentration constant and accounting for variability among studies and species.

Prediction figure

#Overlaying glmm results on figure, but prediction line is jagged...
pred_modDSb11_log <- predict(modDSb11_log, newdata=largeDS2, type="response")
largeDS3 <- cbind(largeDS2, pred_modDSb11_log)

largeDS_plot <- ggplot(data = largeDS3) +
    geom_point(mapping = aes(
      x = Sed_level_stand_mg, 
      y = Binary_large_necroses,
      color = Ref)) +
    labs(x = expression("Sediment exposure concentration (mg/cm"^"2"*"/day)"), 
         y = "Large necroses due to sediment exposure?",
         color = "Study") +
    scale_x_log10(limits=c(1,1000),breaks=c(0.01,0.1,1,10,100,1000),
                  label=c("0.01","0.1","1","10","100","1000")) +
    geom_line(aes(x = Sed_level_stand_mg, y = pred_modDSb11_log), inherit.aes=FALSE) +
    theme_bw()

largeDS_plot

That plot is confusing to interpret. Let’s see if I can plot the average marginal probability, i.e., the average change in probability of the outcome across the range of the predictor of interest. This is described in some detail at the following, useful website: https://stats.idre.ucla.edu/r/dae/mixed-effects-logistic-regression/

By exposure concentration
jvalues <- with(largeDS2, seq(from = min(log10_conc), to = max(log10_conc), length.out = 100))

# calculate predicted probabilities and store in a list
pp <- lapply(jvalues, function(j) {
    largeDS2$log10_conc <- j
    predict(modDSb11_log, newdata = largeDS2, type = "response")
})

# average marginal predicted probability across a few different sediment levels
#sapply(pp[c(1, 20, 40, 60, 80, 100)], mean)
## 
# get the means with lower and upper quartiles
plotdat <- t(sapply(pp, function(x) {
    c(M = mean(x), med = median(x), quantile(x, c(0.25, 0.75)))
}))

# add in log10_conc values and convert to data frame
plotdat <- as.data.frame(cbind(plotdat, jvalues))
plotdat <- plotdat %>% mutate(Sed_level_linear=10^jvalues)

# better names and show the first few rows
colnames(plotdat) <- c("PredictedProbability", "MedianProbability", 
                       "LowerQuantile", "UpperQuantile", "log10_conc", "Sed_conc")
#head(plotdat)

# plot average marginal predicted probabilities
ggplot(plotdat, aes(x = log10_conc)) + 
  geom_line(aes(y = PredictedProbability), size = 2) +
  geom_line(aes(y = MedianProbability), size = 0.5) +
  ylim(c(0, 1))

ggplot(plotdat, aes(x = log10_conc)) + 
  geom_ribbon(aes(ymin = LowerQuantile, ymax = UpperQuantile), alpha = 0.15) +
  geom_line(aes(y = PredictedProbability), size = 2) +
  geom_line(aes(y = MedianProbability), size = 0.5) + 
  ylim(c(0, 1))

ggplot(plotdat, aes(x = Sed_conc, y = PredictedProbability)) + 
  geom_ribbon(aes(ymin = LowerQuantile, ymax = UpperQuantile), alpha = 0.15) +
  geom_line(aes(y = PredictedProbability), size = 2) +
  geom_line(aes(y = MedianProbability), size = 0.5) + 
  ylim(c(0, 1)) +
  scale_x_log10()

#Overlaying average marginal predicted probabilities on figure
#by Ref
largeDS_plot2 <- ggplot() +
    geom_point(data = largeDS2, mapping = aes(
      x = Sed_level_stand_mg, 
      y = Binary_large_necroses)) +
    labs(x = expression("Sediment exposure concentration (mg/cm"^"2"*"/day)"), 
         y = "Predicted probability of large\nnecroses due to sediment exposure",
         linetype = "Predicted\nProbability") +
    scale_x_log10(limits=c(1,1000), breaks=c(0.01,0.1,1,10,100,1000,loaelDS2),
                  label=c("0.01","0.1","1","10","100","1000",round(loaelDS2,digits=1))) +
    geom_ribbon(data = plotdat, aes(x = Sed_conc, y = PredictedProbability, 
                                    ymin = LowerQuantile, ymax = UpperQuantile), 
                alpha = 0.15) +
    geom_line(data = plotdat, aes(x = Sed_conc, y = PredictedProbability, 
                                  linetype = "twodash")) +
    geom_line(data = plotdat, aes(x = Sed_conc, y = MedianProbability, 
                                  linetype = "solid")) +
    geom_vline(xintercept=loaelDS2, linetype="dashed", color = "red") +
    theme_bw() +
    scale_linetype_manual(values=c("twodash", "solid"), labels = c("Median","Mean"))

largeDS_plot2

By exposure duration
jvalues <- with(largeDS2, seq(from = min(log10_dur), to = max(log10_dur), length.out = 100))

# calculate predicted probabilities and store in a list
pp <- lapply(jvalues, function(j) {
    largeDS2$log10_dur <- j
    predict(modDSb11_log, newdata = largeDS2, type = "response")
})

# average marginal predicted probability across a few different sediment levels
#sapply(pp[c(1, 20, 40, 60, 80, 100)], mean)
## 
# get the means with lower and upper quartiles
plotdat <- t(sapply(pp, function(x) {
    c(M = mean(x), med = median(x), quantile(x, c(0.25, 0.75)))
}))

# add in log10_dur values and convert to data frame
plotdat <- as.data.frame(cbind(plotdat, jvalues))
plotdat <- plotdat %>% mutate(Sed_level_linear=10^jvalues)

# better names and show the first few rows
colnames(plotdat) <- c("PredictedProbability", "MedianProbability", 
                       "LowerQuantile", "UpperQuantile", "log10_dur", "Sed_dur")
#head(plotdat)

# plot average marginal predicted probabilities
ggplot(plotdat, aes(x = log10_dur)) + 
  geom_line(aes(y = PredictedProbability), size = 2) +
  geom_line(aes(y = MedianProbability), size = 0.5) +
  ylim(c(0, 1))

ggplot(plotdat, aes(x = log10_dur)) + 
  geom_ribbon(aes(ymin = LowerQuantile, ymax = UpperQuantile), alpha = 0.15) +
  geom_line(aes(y = PredictedProbability), size = 2) +
  geom_line(aes(y = MedianProbability), size = 0.5) + 
  ylim(c(0, 1))

ggplot(plotdat, aes(x = Sed_dur, y = PredictedProbability)) + 
  geom_ribbon(aes(ymin = LowerQuantile, ymax = UpperQuantile), alpha = 0.15) +
  geom_line(aes(y = PredictedProbability), size = 2) +
  geom_line(aes(y = MedianProbability), size = 0.5) + 
  ylim(c(0, 1)) +
  scale_x_log10()

#Overlaying average marginal predicted probabilities on figure
#by Ref
largeDS_plot2 <- ggplot() +
    geom_point(data = largeDS2, mapping = aes(
      x = Sed_exposure_days, 
      y = Binary_large_necroses)) +
    labs(x = expression("Sediment exposure duration (days)"), 
         y = "Predicted probability of large\nnecroses due to sediment exposure",
         linetype = "Predicted\nProbability") +
    scale_x_log10(limits=c(0.1,max(largeDS2$Sed_exposure_days)), breaks=c(0.01,0.1,1,10,100,1000),
                  label=c("0.01","0.1","1","10","100","1000")) +
    geom_ribbon(data = plotdat, aes(x = Sed_dur, y = PredictedProbability, 
                                    ymin = LowerQuantile, ymax = UpperQuantile), 
                alpha = 0.15) +
    geom_line(data = plotdat, aes(x = Sed_dur, y = PredictedProbability, 
                                  linetype = "twodash")) +
    geom_line(data = plotdat, aes(x = Sed_dur, y = MedianProbability, 
                                  linetype = "solid")) +
    theme_bw() +
    scale_linetype_manual(values=c("twodash", "solid"), labels = c("Median","Mean"))

largeDS_plot2

Total mortality

##     n
## 1 509
## Data: totalDS2
## Models:
## modDSc6: Binary_death_total ~ Sed_level_stand_mg * Sed_exposure_days + 
## modDSc6:     (1 | Gsp)
## modDSc9: Binary_death_total ~ Sed_level_stand_mg * Sed_exposure_days + 
## modDSc9:     (1 | Ref)
## modDSc12: Binary_death_total ~ Sed_level_stand_mg * Sed_exposure_days + 
## modDSc12:     (1 | Gsp) + (1 | Ref)
##          npar    AIC    BIC  logLik deviance  Chisq Df Pr(>Chisq)    
## modDSc6     5 306.66 327.82 -148.33   296.66                         
## modDSc9     5 392.75 413.91 -191.37   382.75  0.000  0          1    
## modDSc12    6 295.74 321.13 -141.87   283.74 99.009  1     <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## Data: totalDS2
## Models:
## modDSc10: Binary_death_total ~ Sed_level_stand_mg + (1 | Gsp) + (1 | Ref)
## modDSc11: Binary_death_total ~ Sed_level_stand_mg + Sed_exposure_days + 
## modDSc11:     (1 | Gsp) + (1 | Ref)
## modDSc12: Binary_death_total ~ Sed_level_stand_mg * Sed_exposure_days + 
## modDSc12:     (1 | Gsp) + (1 | Ref)
##          npar    AIC    BIC  logLik deviance  Chisq Df Pr(>Chisq)    
## modDSc10    4 308.19 325.12 -150.09   300.19                         
## modDSc11    5 295.52 316.69 -142.76   285.52 14.664  1  0.0001285 ***
## modDSc12    6 295.74 321.13 -141.87   283.74  1.784  1  0.1816564    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## Data: totalDS2
## Models:
## modDSc4: Binary_death_total ~ Sed_level_stand_mg + (1 | Gsp)
## modDSc7: Binary_death_total ~ Sed_level_stand_mg + (1 | Ref)
## modDSc10: Binary_death_total ~ Sed_level_stand_mg + (1 | Gsp) + (1 | Ref)
##          npar    AIC    BIC  logLik deviance  Chisq Df Pr(>Chisq)    
## modDSc4     3 319.21 331.90 -156.60   313.21                         
## modDSc7     3 396.32 409.02 -195.16   390.32  0.000  0          1    
## modDSc10    4 308.19 325.12 -150.09   300.19 90.131  1     <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: binomial  ( logit )
## Formula: Binary_death_total ~ log10_conc + log10_dur + (1 | Gsp) + (1 |  
##     Ref)
##    Data: totalDS2
## 
##      AIC      BIC   logLik deviance df.resid 
##    255.0    276.2   -122.5    245.0      504 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.2021 -0.0137 -0.0021  0.0000  7.6530 
## 
## Random effects:
##  Groups Name        Variance Std.Dev.
##  Gsp    (Intercept) 280.27   16.741  
##  Ref    (Intercept)  42.47    6.517  
## Number of obs: 509, groups:  Gsp, 75; Ref, 23
## 
## Fixed effects:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  -38.909      6.755  -5.760 8.42e-09 ***
## log10_conc     6.103      1.500   4.069 4.73e-05 ***
## log10_dur      7.855      1.711   4.590 4.44e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Effect estimates

## [1] 0.9685658
##    
## p     0   1
##   0 419   9
##   1   7  74
## Setting levels: control = 0, case = 1
## Setting direction: controls < cases
## Area under the curve: 0.9905
##                      R2m        R2c
## theoretical 0.0683368385 0.99059879
## delta       0.0009323897 0.01351576

## $Gsp

## 
## $Ref

##  Groups Name        Std.Dev.
##  Gsp    (Intercept) 16.741  
##  Ref    (Intercept)  6.517
## [1] 676.5725

##                      Est           LL           UL
## (Intercept) 1.264562e-17 2.247393e-23 7.115428e-12
## log10_conc  4.471818e+02 2.363854e+01 8.459556e+03
## log10_dur   2.577456e+03 9.003714e+01 7.378375e+04

For every 10-fold increase in exposure concentration of deposited sediment, the odds of total colony mortality increase by 447 times (95% CI 24, 8459; GLMM z = 4.069, p < 0.0001), after accounting for exposure duration and variability among studies and species. For every 10-fold change in exposure duration, the odds of total colony mortality increase by 2577 times (95% CI 90, 73,784; GLMM z = 4.590, p = 0.00001), after accounting for exposure concentration and variability among studies and species.

Prediction figure

#Overlaying glmm results on figure, but prediction line is jagged...
pred_modDSc11_log <- predict(modDSc11_log, newdata=totalDS2, type="response")
totalDS3 <- cbind(totalDS2, pred_modDSc11_log)

totalDS_plot <- ggplot(data = totalDS3) +
    geom_point(mapping = aes(
      x = Sed_level_stand_mg, 
      y = Binary_death_total,
      color = Ref)) +
    labs(x = expression("Sediment exposure concentration (mg/cm"^"2"*"/day)"), 
         y = "Total mortality due to sediment exposure?",
         color = "Study") +
    scale_x_log10(limits=c(0.01,1000),breaks=c(0.01,0.1,1,10,100,1000),
                  label=c("0.01","0.1","1","10","100","1000")) +
    geom_line(aes(x = Sed_level_stand_mg, y = pred_modDSc11_log), inherit.aes=FALSE) +
    theme_bw()

totalDS_plot

That plot is confusing to interpret. Let’s see if I can plot the average marginal probability, i.e., the average change in probability of the outcome across the range of the predictor of interest. This is described in some detail at the following, useful website: https://stats.idre.ucla.edu/r/dae/mixed-effects-logistic-regression/

By exposure concentration
jvalues <- with(totalDS2, seq(from = min(log10_conc), to = max(log10_conc), length.out = 100))

# calculate predicted probabilities and store in a list
pp <- lapply(jvalues, function(j) {
    totalDS2$log10_conc <- j
    predict(modDSc11_log, newdata = totalDS2, type = "response")
})

# average marginal predicted probability across a few different sediment levels
#sapply(pp[c(1, 20, 40, 60, 80, 100)], mean)
## 
# get the means with lower and upper quartiles
plotdat <- t(sapply(pp, function(x) {
    c(M = mean(x), med = median(x), quantile(x, c(0.25, 0.75)))
}))

# add in log10_conc values and convert to data frame
plotdat <- as.data.frame(cbind(plotdat, jvalues))
plotdat <- plotdat %>% mutate(Sed_level_linear=10^jvalues)

# better names and show the first few rows
colnames(plotdat) <- c("PredictedProbability", "MedianProbability", 
                       "LowerQuantile", "UpperQuantile", "log10_conc", "Sed_conc")
#head(plotdat)

# plot average marginal predicted probabilities
ggplot(plotdat, aes(x = log10_conc)) + 
  geom_line(aes(y = PredictedProbability), size = 2) +
  geom_line(aes(y = MedianProbability), size = 0.5) +
  ylim(c(0, 1))

ggplot(plotdat, aes(x = log10_conc)) + 
  geom_ribbon(aes(ymin = LowerQuantile, ymax = UpperQuantile), alpha = 0.15) +
  geom_line(aes(y = PredictedProbability), size = 2) +
  geom_line(aes(y = MedianProbability), size = 0.5) + 
  ylim(c(0, 1))

ggplot(plotdat, aes(x = Sed_conc, y = PredictedProbability)) + 
  geom_ribbon(aes(ymin = LowerQuantile, ymax = UpperQuantile), alpha = 0.15) +
  geom_line(aes(y = PredictedProbability), size = 2) +
  geom_line(aes(y = MedianProbability), size = 0.5) + 
  ylim(c(0, 1)) +
  scale_x_log10()

#Overlaying average marginal predicted probabilities on figure
#by Ref
totalDS_plot2 <- ggplot() +
    geom_point(data = totalDS2, mapping = aes(
      x = Sed_level_stand_mg, 
      y = Binary_death_total,
      color = Ref)) +
    labs(x = expression("Sediment exposure concentration (mg/cm"^"2"*"/day)"), 
         y = "Predicted probability of total colony\nmortality due to sediment exposure",
         color = "Study",
         linetype = "Predicted Probability") +
    scale_x_log10(limits=c(0.01,1000), breaks=c(0.01,0.1,1,10,100,1000,loaelDS3),
                  label=c("0.01","0.1","1","10","100","1000",round(loaelDS3,digits=1))) +
    geom_ribbon(data = plotdat, aes(x = Sed_conc, y = PredictedProbability, 
                                    ymin = LowerQuantile, ymax = UpperQuantile), 
                alpha = 0.15) +
    geom_line(data = plotdat, aes(x = Sed_conc, y = PredictedProbability, 
                                  linetype = "twodash")) +
    geom_line(data = plotdat, aes(x = Sed_conc, y = MedianProbability, 
                                  linetype = "solid")) +
    geom_vline(xintercept=loaelDS3, linetype="dashed", color = "red") +
    theme_bw() +
    scale_linetype_manual(values=c("twodash", "solid"), labels = c("Median","Mean"))

totalDS_plot2

By exposure duration
jvalues <- with(totalDS2, seq(from = min(log10_dur), to = max(log10_dur), length.out = 100))

# calculate predicted probabilities and store in a list
pp <- lapply(jvalues, function(j) {
    totalDS2$log10_dur <- j
    predict(modDSc11_log, newdata = totalDS2, type = "response")
})

# get the means with lower and upper quartiles
plotdat <- t(sapply(pp, function(x) {
    c(M = mean(x), med = median(x), quantile(x, c(0.25, 0.75)))
}))

# add in log10_dur values and convert to data frame
plotdat <- as.data.frame(cbind(plotdat, jvalues))
plotdat <- plotdat %>% mutate(Sed_level_linear=10^jvalues)

# better names and show the first few rows
colnames(plotdat) <- c("PredictedProbability", "MedianProbability", 
                       "LowerQuantile", "UpperQuantile", "log10_dur", "Sed_dur")
#head(plotdat)

# plot average marginal predicted probabilities
ggplot(plotdat, aes(x = log10_dur)) + 
  geom_line(aes(y = PredictedProbability), size = 2) +
  geom_line(aes(y = MedianProbability), size = 0.5) +
  ylim(c(0, 1))

ggplot(plotdat, aes(x = log10_dur)) + 
  geom_ribbon(aes(ymin = LowerQuantile, ymax = UpperQuantile), alpha = 0.15) +
  geom_line(aes(y = PredictedProbability), size = 2) +
  geom_line(aes(y = MedianProbability), size = 0.5) + 
  ylim(c(0, 1))

ggplot(plotdat, aes(x = Sed_dur, y = PredictedProbability)) + 
  geom_ribbon(aes(ymin = LowerQuantile, ymax = UpperQuantile), alpha = 0.15) +
  geom_line(aes(y = PredictedProbability), size = 2) +
  geom_line(aes(y = MedianProbability), size = 0.5) + 
  ylim(c(0, 1)) +
  scale_x_log10()

#Overlaying average marginal predicted probabilities on figure
#by Ref
totalDS_plot2 <- ggplot() +
    geom_point(data = totalDS2, mapping = aes(
      x = Sed_exposure_days, 
      y = Binary_death_total)) +
    labs(x = expression("Sediment exposure duration (days)"), 
         y = "Predicted probability of total colony\nmortality due to sediment exposure",
         linetype = "Predicted\nProbability") +
    scale_x_log10(limits=c(0.1,max(totalDS2$Sed_exposure_days)), breaks=c(0.01,0.1,1,10,100,1000),
                  label=c("0.01","0.1","1","10","100","1000")) +
    geom_ribbon(data = plotdat, aes(x = Sed_dur, y = PredictedProbability, 
                                    ymin = LowerQuantile, ymax = UpperQuantile), 
                alpha = 0.15) +
    geom_line(data = plotdat, aes(x = Sed_dur, y = PredictedProbability, 
                                  linetype = "twodash")) +
    geom_line(data = plotdat, aes(x = Sed_dur, y = MedianProbability, 
                                  linetype = "solid")) +
    theme_bw() +
    scale_linetype_manual(values=c("twodash", "solid"), labels = c("Median","Mean"))

totalDS_plot2

SUSPENDED SEDIMENT

Small Necroses

##     n
## 1 147
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## Data: smallSS2
## Models:
## modSS6: Binary_small_necroses ~ Sed_level_stand_mg * Sed_exposure_days + 
## modSS6:     (1 | Gsp)
## modSS9: Binary_small_necroses ~ Sed_level_stand_mg * Sed_exposure_days + 
## modSS9:     (1 | Ref)
## modSS12: Binary_small_necroses ~ Sed_level_stand_mg * Sed_exposure_days + 
## modSS12:     (1 | Gsp) + (1 | Ref)
## modSS15: Binary_small_necroses ~ Sed_level_stand_mg * Sed_exposure_days + 
## modSS15:     (1 | Ref_name/Ref)
## modSS18: Binary_small_necroses ~ Sed_level_stand_mg * Sed_exposure_days + 
## modSS18:     (1 | Gsp) + (1 | Ref_name)
## modSS21: Binary_small_necroses ~ Sed_level_stand_mg * Sed_exposure_days + 
## modSS21:     (1 | Gsp) + (1 | Ref_name/Ref)
##         npar    AIC    BIC  logLik deviance  Chisq Df Pr(>Chisq)    
## modSS6     5 108.19 123.14 -49.094   98.189                         
## modSS9     5 126.29 141.24 -58.145  116.290  0.000  0          1    
## modSS12    6 110.19 128.13 -49.094   98.189 18.102  1  2.094e-05 ***
## modSS15    6 128.29 146.23 -58.145  116.290  0.000  0          1    
## modSS18    6 110.19 128.13 -49.094   98.189 18.102  0  < 2.2e-16 ***
## modSS21    7 112.19 133.12 -49.094   98.189  0.000  1          1    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## Data: smallSS2
## Models:
## modSS6: Binary_small_necroses ~ Sed_level_stand_mg * Sed_exposure_days + 
## modSS6:     (1 | Gsp)
## modSS12: Binary_small_necroses ~ Sed_level_stand_mg * Sed_exposure_days + 
## modSS12:     (1 | Gsp) + (1 | Ref)
## modSS18: Binary_small_necroses ~ Sed_level_stand_mg * Sed_exposure_days + 
## modSS18:     (1 | Gsp) + (1 | Ref_name)
##         npar    AIC    BIC  logLik deviance Chisq Df Pr(>Chisq)
## modSS6     5 108.19 123.14 -49.094   98.189                    
## modSS12    6 110.19 128.13 -49.094   98.189     0  1          1
## modSS18    6 110.19 128.13 -49.094   98.189     0  0          1
## Data: smallSS2
## Models:
## modSS4: Binary_small_necroses ~ Sed_level_stand_mg + (1 | Gsp)
## modSS5: Binary_small_necroses ~ Sed_level_stand_mg + Sed_exposure_days + 
## modSS5:     (1 | Gsp)
## modSS6: Binary_small_necroses ~ Sed_level_stand_mg * Sed_exposure_days + 
## modSS6:     (1 | Gsp)
##        npar    AIC    BIC  logLik deviance   Chisq Df Pr(>Chisq)    
## modSS4    3 135.52 144.49 -64.759  129.518                          
## modSS5    4 106.75 118.71 -49.375   98.750 30.7680  1  2.908e-08 ***
## modSS6    5 108.19 123.14 -49.094   98.189  0.5607  1      0.454    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## Data: smallSS2
## Models:
## modSS10: Binary_small_necroses ~ Sed_level_stand_mg + (1 | Gsp) + (1 | 
## modSS10:     Ref)
## modSS11: Binary_small_necroses ~ Sed_level_stand_mg + Sed_exposure_days + 
## modSS11:     (1 | Gsp) + (1 | Ref)
## modSS12: Binary_small_necroses ~ Sed_level_stand_mg * Sed_exposure_days + 
## modSS12:     (1 | Gsp) + (1 | Ref)
##         npar    AIC    BIC  logLik deviance   Chisq Df Pr(>Chisq)    
## modSS10    4 137.45 149.41 -64.724  129.447                          
## modSS11    5 108.75 123.70 -49.375   98.750 30.6976  1  3.015e-08 ***
## modSS12    6 110.19 128.13 -49.094   98.189  0.5607  1      0.454    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## Data: smallSS2
## Models:
## modSS16: Binary_small_necroses ~ Sed_level_stand_mg + (1 | Gsp) + (1 | 
## modSS16:     Ref_name)
## modSS17: Binary_small_necroses ~ Sed_level_stand_mg + Sed_exposure_days + 
## modSS17:     (1 | Gsp) + (1 | Ref_name)
## modSS18: Binary_small_necroses ~ Sed_level_stand_mg * Sed_exposure_days + 
## modSS18:     (1 | Gsp) + (1 | Ref_name)
##         npar    AIC    BIC  logLik deviance   Chisq Df Pr(>Chisq)    
## modSS16    4 137.45 149.41 -64.724  129.447                          
## modSS17    5 108.75 123.70 -49.375   98.750 30.6976  1  3.015e-08 ***
## modSS18    6 110.19 128.13 -49.094   98.189  0.5607  1      0.454    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## Data: smallSS2
## Models:
## modSS5: Binary_small_necroses ~ Sed_level_stand_mg + Sed_exposure_days + 
## modSS5:     (1 | Gsp)
## modSS11: Binary_small_necroses ~ Sed_level_stand_mg + Sed_exposure_days + 
## modSS11:     (1 | Gsp) + (1 | Ref)
## modSS17: Binary_small_necroses ~ Sed_level_stand_mg + Sed_exposure_days + 
## modSS17:     (1 | Gsp) + (1 | Ref_name)
##         npar    AIC    BIC  logLik deviance Chisq Df Pr(>Chisq)
## modSS5     4 106.75 118.71 -49.375    98.75                    
## modSS11    5 108.75 123.70 -49.375    98.75     0  1          1
## modSS17    5 108.75 123.70 -49.375    98.75     0  0          1
## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: binomial  ( logit )
## Formula: Binary_small_necroses ~ Sed_level_stand_mg + Sed_exposure_days +  
##     (1 | Gsp)
##    Data: smallSS2
## 
##      AIC      BIC   logLik deviance df.resid 
##    106.7    118.7    -49.4     98.7      143 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -1.9893 -0.3421 -0.1326 -0.0291  3.8738 
## 
## Random effects:
##  Groups Name        Variance Std.Dev.
##  Gsp    (Intercept) 3.44     1.855   
## Number of obs: 147, groups:  Gsp, 8
## 
## Fixed effects:
##                    Estimate Std. Error z value Pr(>|z|)    
## (Intercept)        -8.06121    1.75339  -4.598 4.28e-06 ***
## Sed_level_stand_mg  0.03325    0.00998   3.331 0.000864 ***
## Sed_exposure_days   0.07866    0.01806   4.355 1.33e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Effect estimates

## [1] 0.8707483
##    
## p     0   1
##   0 108  11
##   1   8  20
## Setting levels: control = 0, case = 1
## Setting direction: controls < cases
## Area under the curve: 0.9189
##                   R2m       R2c
## theoretical 0.3536090 0.6839914
## delta       0.2701203 0.5224978

## $Gsp

##  Groups Name        Std.Dev.
##  Gsp    (Intercept) 1.8546
## [1] 6.389086

##                            Est           LL         UL
## (Intercept)        0.003743975 0.0003457821 0.04053811
## Sed_level_stand_mg 1.023311817 1.0095310571 1.03728069
## Sed_exposure_days  1.056033682 1.0304395187 1.08226356
##                             Est           LL           UL
## (Intercept)        8.685337e-09 3.177978e-12 2.373682e-05
## Sed_level_stand_mg 1.079558e+00 1.032013e+00 1.129292e+00
## Sed_exposure_days  1.198549e+00 1.104739e+00 1.300325e+00

For every 10-fold increase in exposure concentration of suspended sediment, the odds of developing small tissue necroses increase by 1.08 times (95% CI 1.03, 1.13; GLMM z = 3.331, p = 0.0008), while holding exposure duration constant and after accounting for variability among species. For every 10-fold increase in exposure duration, the odds of small tissue necroses increase by 1.20 times (95% CI 1.10, 1.30; GLMM z = 4.355, p < 0.0001), while holding exposure concentration constant and after accounting for variability among species.

Prediction figure

#Overlaying glmm results on figure, but prediction line is jagged...
pred_modSS5 <- predict(modSS5, newdata=smallSS2, type="response")
smallSS3 <- cbind(smallSS2, pred_modSS5)

smallSS_plot <- ggplot(data = smallSS3) +
    geom_point(mapping = aes(
      x = Sed_level_stand_mg, 
      y = Binary_small_necroses,
      color = Gsp)) +
    labs(x = "Sediment exposure concentration (mg/L)", 
         y = "Small necroses due to sediment exposure?",
         color = "Species") +
    scale_x_log10(limits=c(1,max(smallSS2$Sed_level_stand_mg)),breaks=c(0.01,0.1,1,10,100,1000),
                  label=c("0.01","0.1","1","10","100","1000")) +
    geom_line(aes(x = Sed_level_stand_mg, y = pred_modSS5), inherit.aes=FALSE) +
    theme_bw()

smallSS_plot

That plot is confusing to interpret. Let’s see if I can plot the average marginal probability, i.e., the average change in probability of the outcome across the range of the predictor of interest. This is described in some detail at the following, useful website: https://stats.idre.ucla.edu/r/dae/mixed-effects-logistic-regression/

By sediment exposure concentration
jvalues <- with(smallSS2, seq(from = min(Sed_level_stand_mg), to = max(Sed_level_stand_mg), length.out = 100))

# calculate predicted probabilities and store in a list
pp <- lapply(jvalues, function(j) {
    smallSS2$Sed_level_stand_mg <- j
    predict(modSS5, newdata = smallSS2, type = "response")
})

# get the means with lower and upper quartiles
plotdat <- t(sapply(pp, function(x) {
    c(M = mean(x), med = median(x), quantile(x, c(0.25, 0.75)))
}))

# add in Sed_level_stand_mg values and convert to data frame
plotdat <- as.data.frame(cbind(plotdat, jvalues))

# better names and show the first few rows
colnames(plotdat) <- c("PredictedProbability", "MedianProbability", 
                       "LowerQuantile", "UpperQuantile", "Sed_conc")
#head(plotdat)

# plot average marginal predicted probabilities
ggplot(plotdat, aes(x = Sed_conc)) + 
  geom_line(aes(y = PredictedProbability), size = 2) +
  geom_line(aes(y = MedianProbability), size = 0.5) +
  ylim(c(0, 1))

ggplot(plotdat, aes(x = Sed_conc)) + 
  geom_ribbon(aes(ymin = LowerQuantile, ymax = UpperQuantile), alpha = 0.15) +
  geom_line(aes(y = PredictedProbability), size = 2) +
  geom_line(aes(y = MedianProbability), size = 0.5) + 
  ylim(c(0, 1))

#Overlaying average marginal predicted probabilities on figure
#by Ref
smallSS_plot2 <- ggplot() +
    geom_point(data = smallSS2, mapping = aes(
      x = Sed_level_stand_mg, 
      y = Binary_small_necroses)) +
    labs(x = "Sediment exposure concentration (mg/L)", 
         y = "Predicted probability of small\nnecroses due to sediment exposure",
         linetype = "Predicted\nProbability") +
    scale_x_continuous(limits=c(0,125), breaks=c(0,25,50,75,100,125,loaelSS1), 
                       labels=c("0","25","50","75","100","125",round(loaelSS1,digits=1))) +
    geom_ribbon(data = plotdat, aes(x = Sed_conc, y = PredictedProbability, 
                                    ymin = LowerQuantile, ymax = UpperQuantile), 
                alpha = 0.15) +
    geom_line(data = plotdat, aes(x = Sed_conc, y = PredictedProbability, 
                                  linetype = "twodash")) +
    geom_line(data = plotdat, aes(x = Sed_conc, y = MedianProbability, 
                                  linetype = "solid")) +
    geom_vline(xintercept=loaelSS1, linetype="dashed", color = "red") +
    theme_bw() +
    scale_linetype_manual(values=c("twodash", "solid"), labels = c("Median","Mean"))

smallSS_plot2

By sediment exposure duration
jvalues <- with(smallSS2, seq(from = min(Sed_exposure_days), to = max(Sed_exposure_days), length.out = 100))

# calculate predicted probabilities and store in a list
pp <- lapply(jvalues, function(j) {
    smallSS2$Sed_exposure_days <- j
    predict(modSS5, newdata = smallSS2, type = "response")
})

# get the means with lower and upper quartiles
plotdat <- t(sapply(pp, function(x) {
    c(M = mean(x), med = median(x), quantile(x, c(0.25, 0.75)))
}))

# add in Sed_exposure_days values and convert to data frame
plotdat <- as.data.frame(cbind(plotdat, jvalues))

# better names and show the first few rows
colnames(plotdat) <- c("PredictedProbability", "MedianProbability", 
                       "LowerQuantile", "UpperQuantile", "Sed_dur")
#head(plotdat)

# plot average marginal predicted probabilities
ggplot(plotdat, aes(x = Sed_dur)) + 
  geom_line(aes(y = PredictedProbability), size = 2) +
  geom_line(aes(y = MedianProbability), size = 0.5) +
  ylim(c(0, 1))

ggplot(plotdat, aes(x = Sed_dur)) + 
  geom_ribbon(aes(ymin = LowerQuantile, ymax = UpperQuantile), alpha = 0.15) +
  geom_line(aes(y = PredictedProbability), size = 2) +
  geom_line(aes(y = MedianProbability), size = 0.5) + 
  ylim(c(0, 1))

#Overlaying average marginal predicted probabilities on figure
#by Ref
smallSS_plot3 <- ggplot() +
    geom_point(data = smallSS2, mapping = aes(
      x = Sed_exposure_days, 
      y = Binary_small_necroses,
      color = Gsp)) +
    labs(x = "Sediment exposure duration (days)", 
         y = "Predicted probability of small\nnecroses to sediment exposure",
         color = "Species",
         linetype = "Predicted\nProbability") +
    scale_x_continuous(limits=c(0,max(smallSS2$Sed_exposure_days))) +
    geom_ribbon(data = plotdat, aes(x = Sed_dur, y = PredictedProbability, 
                                    ymin = LowerQuantile, ymax = UpperQuantile), 
                alpha = 0.15) +
    geom_line(data = plotdat, aes(x = Sed_dur, y = PredictedProbability, 
                                  linetype = "twodash")) +
    geom_line(data = plotdat, aes(x = Sed_dur, y = MedianProbability, 
                                  linetype = "solid")) +
    theme_bw() +
    scale_linetype_manual(values=c("twodash", "solid"), labels = c("Median","Mean"))

smallSS_plot3

Large Necroses

##     n
## 1 147
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
##         df      AIC
## modSSb1  2 52.19214
## modSSb2  3 11.54518
## modSSb3  4 13.54518
## Data: largeSS2
## Models:
## modSSb7: Binary_large_necroses ~ Sed_level_stand_mg + (1 | Ref)
## modSSb8: Binary_large_necroses ~ Sed_level_stand_mg + Sed_exposure_days + 
## modSSb8:     (1 | Ref)
## modSSb9: Binary_large_necroses ~ Sed_level_stand_mg * Sed_exposure_days + 
## modSSb9:     (1 | Ref)
##         npar    AIC    BIC   logLik deviance  Chisq Df Pr(>Chisq)    
## modSSb7    3 31.918 40.889 -12.9588  25.9176                         
## modSSb8    4 13.545 25.507  -2.7726   5.5452 20.372  1  6.374e-06 ***
## modSSb9    5 15.545 30.497  -2.7726   5.5452  0.000  1          1    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## Data: largeSS2
## Models:
## modSSb13: Binary_large_necroses ~ Sed_level_stand_mg + (1 | Ref_name/Ref)
## modSSb14: Binary_large_necroses ~ Sed_level_stand_mg + Sed_exposure_days + 
## modSSb14:     (1 | Ref_name/Ref)
## modSSb15: Binary_large_necroses ~ Sed_level_stand_mg * Sed_exposure_days + 
## modSSb15:     (1 | Ref_name/Ref)
##          npar    AIC    BIC   logLik deviance  Chisq Df Pr(>Chisq)    
## modSSb13    4 33.918 45.879 -12.9588  25.9176                         
## modSSb14    5 15.545 30.497  -2.7726   5.5452 20.372  1  6.374e-06 ***
## modSSb15    6 17.545 35.488  -2.7726   5.5452  0.000  1          1    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## Data: largeSS2
## Models:
## modSSb4: Binary_large_necroses ~ Sed_level_stand_mg + (1 | Gsp)
## modSSb8: Binary_large_necroses ~ Sed_level_stand_mg + Sed_exposure_days + 
## modSSb8:     (1 | Ref)
## modSSb10: Binary_large_necroses ~ Sed_level_stand_mg + (1 | Gsp) + (1 | 
## modSSb10:     Ref)
## modSSb16: Binary_large_necroses ~ Sed_level_stand_mg + (1 | Gsp) + (1 | 
## modSSb16:     Ref_name)
## modSSb14: Binary_large_necroses ~ Sed_level_stand_mg + Sed_exposure_days + 
## modSSb14:     (1 | Ref_name/Ref)
## modSSb19: Binary_large_necroses ~ Sed_level_stand_mg + (1 | Gsp) + (1 | 
## modSSb19:     Ref_name/Ref)
##          npar    AIC    BIC   logLik deviance  Chisq Df Pr(>Chisq)    
## modSSb4     3 36.784 45.756 -15.3921  30.7842                         
## modSSb8     4 13.545 25.507  -2.7726   5.5452 25.239  1  5.065e-07 ***
## modSSb10    4 33.803 45.765 -12.9014  25.8028  0.000  0          1    
## modSSb16    4 33.803 45.765 -12.9014  25.8028  0.000  0          1    
## modSSb14    5 15.545 30.497  -2.7726   5.5452 20.258  1  6.768e-06 ***
## modSSb19    5 35.803 50.755 -12.9014  25.8028  0.000  0          1    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## Data: largeSS2
## Models:
## modSSb8: Binary_large_necroses ~ Sed_level_stand_mg + Sed_exposure_days + 
## modSSb8:     (1 | Ref)
## modSSb14: Binary_large_necroses ~ Sed_level_stand_mg + Sed_exposure_days + 
## modSSb14:     (1 | Ref_name/Ref)
##          npar    AIC    BIC  logLik deviance Chisq Df Pr(>Chisq)
## modSSb8     4 13.545 25.507 -2.7726   5.5452                    
## modSSb14    5 15.545 30.497 -2.7726   5.5452     0  1          1
## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: binomial  ( logit )
## Formula: Binary_large_necroses ~ Sed_level_stand_mg + Sed_exposure_days +  
##     (1 | Ref)
##    Data: largeSS2
## 
##      AIC      BIC   logLik deviance df.resid 
##     13.5     25.5     -2.8      5.5      143 
## 
## Scaled residuals: 
##    Min     1Q Median     3Q    Max 
##     -1      0      0      0      1 
## 
## Random effects:
##  Groups Name        Variance Std.Dev.
##  Ref    (Intercept) 0        0       
## Number of obs: 147, groups:  Ref, 4
## 
## Fixed effects:
##                    Estimate Std. Error z value Pr(>|z|)  
## (Intercept)        -779.798    968.212  -0.805   0.4206  
## Sed_level_stand_mg    2.026     31.654   0.064   0.9490  
## Sed_exposure_days     8.582      3.984   2.154   0.0312 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## convergence code: 0
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: binomial  ( logit )
## Formula: Binary_large_necroses ~ log10_conc + log10_dur + (1 | Ref)
##    Data: largeSS2
## 
##      AIC      BIC   logLik deviance df.resid 
##     13.5     25.5     -2.8      5.5      143 
## 
## Scaled residuals: 
##    Min     1Q Median     3Q    Max 
##     -1      0      0      0      1 
## 
## Random effects:
##  Groups Name        Variance  Std.Dev.
##  Ref    (Intercept) 3.193e-17 5.65e-09
## Number of obs: 147, groups:  Ref, 4
## 
## Fixed effects:
##               Estimate Std. Error z value Pr(>|z|)
## (Intercept)    -1347.0 49256351.2       0        1
## log10_conc        85.8 16165394.5       0        1
## log10_dur        634.7 15817710.9       0        1
## convergence code: 0
## boundary (singular) fit: see ?isSingular
## 
## Call:
## glm(formula = Binary_large_necroses ~ Sed_level_stand_mg + Sed_exposure_days, 
##     family = binomial(link = "logit"), data = largeSS2)
## 
## Deviance Residuals: 
##    Min      1Q  Median      3Q     Max  
## -1.177   0.000   0.000   0.000   1.177  
## 
## Coefficients:
##                      Estimate Std. Error z value Pr(>|z|)
## (Intercept)          -390.084 100469.392  -0.004    0.997
## Sed_level_stand_mg      1.077    314.538   0.003    0.997
## Sed_exposure_days       4.271   1100.945   0.004    0.997
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 50.1358  on 146  degrees of freedom
## Residual deviance:  5.5452  on 144  degrees of freedom
## AIC: 11.545
## 
## Number of Fisher Scoring iterations: 25

Effect estimates

## [1] 0.9863946
##    
## p     0   1
##   0 141   2
##   1   0   4
## Setting levels: control = 0, case = 1
## Setting direction: controls < cases
## Area under the curve: 0.9976
##                   R2m       R2c
## theoretical 0.9996563 0.9996563
## delta       0.9973381 0.9973381

There is no significant relationship between sediment exposure and the odds of developing large tissue necroses (concentration: GLMM z = 0.003, p = 0.997; duration: GLMM z = 0.004, p = 0.997).

Prediction figure

#Overlaying glmm results on figure, but prediction line is jagged...
pred_modSSb2 <- predict(modSSb2, newdata=largeSS2, type="response")
largeSS3 <- cbind(largeSS2, pred_modSSb2)

largeSS_plot <- ggplot(data = largeSS3) +
    geom_point(mapping = aes(
      x = Sed_level_stand_mg, 
      y = Binary_large_necroses,
      color = Ref)) +
    labs(x = "Sediment exposure concentration (mg/L)", 
         y = "Large necroses due to sediment exposure?",
         color = "Study") +
    scale_x_continuous(limits=c(0,max(largeSS2$Sed_level_stand_mg))) +
    geom_line(aes(x = Sed_level_stand_mg, y = pred_modSSb2), inherit.aes=FALSE) +
    theme_bw()

largeSS_plot

That plot is confusing to interpret. Let’s see if I can plot the average marginal probability, i.e., the average change in probability of the outcome across the range of the predictor of interest. This is described in some detail at the following, useful website: https://stats.idre.ucla.edu/r/dae/mixed-effects-logistic-regression/

By exposure concentration
jvalues <- with(largeSS2, seq(from = min(Sed_level_stand_mg), to = max(Sed_level_stand_mg), length.out = 100))

# calculate predicted probabilities and store in a list
pp <- lapply(jvalues, function(j) {
    largeSS2$Sed_level_stand_mg <- j
    predict(modSSb2, newdata = largeSS2, type = "response")
})

# get the means with lower and upper quartiles
plotdat <- t(sapply(pp, function(x) {
    c(M = mean(x), med = median(x), quantile(x, c(0.25, 0.75)))
}))

# add in Sed_level_stand_mg values and convert to data frame
plotdat <- as.data.frame(cbind(plotdat, jvalues))

# better names and show the first few rows
colnames(plotdat) <- c("PredictedProbability", "MedianProbability", 
                       "LowerQuantile", "UpperQuantile", "Sed_conc")
#head(plotdat)

# plot average marginal predicted probabilities
ggplot(plotdat, aes(x = Sed_conc)) + 
  geom_line(aes(y = PredictedProbability), size = 2) +
  geom_line(aes(y = MedianProbability), size = 0.5) +
  ylim(c(0, 1))

ggplot(plotdat, aes(x = Sed_conc)) + 
  geom_ribbon(aes(ymin = LowerQuantile, ymax = UpperQuantile), alpha = 0.15) +
  geom_line(aes(y = PredictedProbability), size = 2) +
  geom_line(aes(y = MedianProbability), size = 0.5) + 
  ylim(c(0, 1))

#Overlaying average marginal predicted probabilities on figure
#by Ref
largeSS_plot2 <- ggplot() +
    geom_point(data = largeSS2, mapping = aes(
      x = Sed_level_stand_mg, 
      y = Binary_large_necroses,
      color = Ref)) +
    labs(x = "Sediment exposure concentration (mg/L)", 
         y = "Predicted probability of large\nnecroses due to sediment exposure",
         color = "Binary data\nby Study",
         linetype = "Predicted\nProbability") +
    scale_x_continuous(limits=c(0,max(largeSS2$Sed_level_stand_mg)),
                       breaks = c(0,25,50,75,100,125,loaelSS2),
                       labels = c("0","25","50","75","100","125",round(loaelSS2,digits=1))) +
    geom_ribbon(data = plotdat, aes(x = Sed_conc, y = PredictedProbability, 
                                    ymin = LowerQuantile, ymax = UpperQuantile), 
                alpha = 0.15) +
    geom_line(data = plotdat, aes(x = Sed_conc, y = PredictedProbability, 
                                  linetype = "twodash")) +
    geom_line(data = plotdat, aes(x = Sed_conc, y = MedianProbability, 
                                  linetype = "solid")) +
    geom_vline(xintercept=loaelSS2, linetype="dashed", color = "red") +
    theme_bw() +
    scale_linetype_manual(values=c("twodash", "solid"), labels = c("Median","Mean"))

largeSS_plot2

By exposure duration
jvalues <- with(largeSS2, seq(from = min(Sed_exposure_days), to = max(Sed_exposure_days), length.out = 100))

# calculate predicted probabilities and store in a list
pp <- lapply(jvalues, function(j) {
    largeSS2$Sed_exposure_days <- j
    predict(modSSb2, newdata = largeSS2, type = "response")
})

# get the means with lower and upper quartiles
plotdat <- t(sapply(pp, function(x) {
    c(M = mean(x), med = median(x), quantile(x, c(0.25, 0.75)))
}))

# add in Sed_exposure_days values and convert to data frame
plotdat <- as.data.frame(cbind(plotdat, jvalues))

# better names and show the first few rows
colnames(plotdat) <- c("PredictedProbability", "MedianProbability", 
                       "LowerQuantile", "UpperQuantile", "Sed_dur")
#head(plotdat)

# plot average marginal predicted probabilities
ggplot(plotdat, aes(x = Sed_dur)) + 
  geom_line(aes(y = PredictedProbability), size = 2) +
  geom_line(aes(y = MedianProbability), size = 0.5) +
  ylim(c(0, 1))

ggplot(plotdat, aes(x = Sed_dur)) + 
  geom_ribbon(aes(ymin = LowerQuantile, ymax = UpperQuantile), alpha = 0.15) +
  geom_line(aes(y = PredictedProbability), size = 2) +
  geom_line(aes(y = MedianProbability), size = 0.5) + 
  ylim(c(0, 1))

#Overlaying average marginal predicted probabilities on figure
#by Ref
largeSS_plot2 <- ggplot() +
    geom_point(data = largeSS2, mapping = aes(
      x = Sed_exposure_days, 
      y = Binary_large_necroses)) +
    labs(x = expression("Sediment exposure duration (days)"), 
         y = "Predicted probability of large\nnecroses due to sediment exposure",
         linetype = "Predicted\nProbability") +
    geom_ribbon(data = plotdat, aes(x = Sed_dur, y = PredictedProbability, 
                                    ymin = LowerQuantile, ymax = UpperQuantile), 
                alpha = 0.15) +
    geom_line(data = plotdat, aes(x = Sed_dur, y = PredictedProbability, 
                                  linetype = "twodash")) +
    geom_line(data = plotdat, aes(x = Sed_dur, y = MedianProbability, 
                                  linetype = "solid")) +
    theme_bw() +
    scale_linetype_manual(values=c("twodash", "solid"), labels = c("Median","Mean"))

largeSS_plot2

Total mortality

##     n
## 1 176
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
##         df      AIC
## modSSc1  2 56.11582
## modSSc2  3 11.54518
## modSSc3  4 13.54518
## Data: totalSS2
## Models:
## modSSc7: Binary_death_total ~ Sed_level_stand_mg + (1 | Ref)
## modSSc8: Binary_death_total ~ Sed_level_stand_mg + Sed_exposure_days + 
## modSSc8:     (1 | Ref)
## modSSc9: Binary_death_total ~ Sed_level_stand_mg * Sed_exposure_days + 
## modSSc9:     (1 | Ref)
##         npar    AIC    BIC   logLik deviance  Chisq Df Pr(>Chisq)    
## modSSc7    3 32.442 41.953 -13.2208  26.4416                         
## modSSc8    4 13.545 26.227  -2.7726   5.5452 20.896  1  4.848e-06 ***
## modSSc9    5 15.545 31.398  -2.7726   5.5452  0.000  1          1    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## Data: totalSS2
## Models:
## modSSc13: Binary_death_total ~ Sed_level_stand_mg + (1 | Ref_name/Ref)
## modSSc15: Binary_death_total ~ Sed_level_stand_mg * Sed_exposure_days + 
## modSSc15:     (1 | Ref_name/Ref)
##          npar    AIC    BIC   logLik deviance  Chisq Df Pr(>Chisq)    
## modSSc13    4 34.441 47.122 -13.2203  26.4406                         
## modSSc15    6 17.545 36.568  -2.7726   5.5452 20.895  2  2.902e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## Data: totalSS2
## Models:
## modSSc4: Binary_death_total ~ Sed_level_stand_mg + (1 | Gsp)
## modSSc8: Binary_death_total ~ Sed_level_stand_mg + Sed_exposure_days + 
## modSSc8:     (1 | Ref)
## modSSc10: Binary_death_total ~ Sed_level_stand_mg + (1 | Gsp) + (1 | Ref)
## modSSc15: Binary_death_total ~ Sed_level_stand_mg * Sed_exposure_days + 
## modSSc15:     (1 | Ref_name/Ref)
##          npar    AIC    BIC   logLik deviance  Chisq Df Pr(>Chisq)    
## modSSc4     3 38.093 47.605 -16.0466   32.093                         
## modSSc8     4 13.545 26.227  -2.7726    5.545 26.548  1  2.571e-07 ***
## modSSc10    4 34.442 47.124 -13.2208   26.442  0.000  0          1    
## modSSc15    6 17.545 36.568  -2.7726    5.545 20.896  2  2.900e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## Data: totalSS2
## Models:
## modSSc8: Binary_death_total ~ Sed_level_stand_mg + Sed_exposure_days + 
## modSSc8:     (1 | Ref)
## modSSc15: Binary_death_total ~ Sed_level_stand_mg * Sed_exposure_days + 
## modSSc15:     (1 | Ref_name/Ref)
##          npar    AIC    BIC  logLik deviance Chisq Df Pr(>Chisq)
## modSSc8     4 13.545 26.227 -2.7726   5.5452                    
## modSSc15    6 17.545 36.568 -2.7726   5.5452     0  2          1
## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: binomial  ( logit )
## Formula: Binary_death_total ~ Sed_level_stand_mg + Sed_exposure_days +  
##     (1 | Ref)
##    Data: totalSS2
## 
##      AIC      BIC   logLik deviance df.resid 
##     13.5     26.2     -2.8      5.5      172 
## 
## Scaled residuals: 
##    Min     1Q Median     3Q    Max 
##     -1      0      0      0      1 
## 
## Random effects:
##  Groups Name        Variance  Std.Dev. 
##  Ref    (Intercept) 5.443e-17 7.378e-09
## Number of obs: 176, groups:  Ref, 8
## 
## Fixed effects:
##                    Estimate Std. Error z value Pr(>|z|)  
## (Intercept)        -862.270    967.393  -0.891   0.3728  
## Sed_level_stand_mg    2.029     31.628   0.064   0.9489  
## Sed_exposure_days     9.562      3.984   2.400   0.0164 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## convergence code: 0
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: binomial  ( logit )
## Formula: Binary_death_total ~ log10_conc + log10_dur + (1 | Ref)
##    Data: totalSS2
## 
##      AIC      BIC   logLik deviance df.resid 
##     13.5     26.2     -2.8      5.5      172 
## 
## Scaled residuals: 
##    Min     1Q Median     3Q    Max 
##     -1      0      0      0      1 
## 
## Random effects:
##  Groups Name        Variance Std.Dev.
##  Ref    (Intercept) 0        0       
## Number of obs: 176, groups:  Ref, 8
## 
## Fixed effects:
##             Estimate Std. Error z value Pr(>|z|)
## (Intercept)  -2676.5     2071.5  -1.292    0.196
## log10_conc      85.2     4512.3   0.019    0.985
## log10_dur     1326.1     3518.0   0.377    0.706
## convergence code: 0
## boundary (singular) fit: see ?isSingular
## 
## Call:
## glm(formula = Binary_death_total ~ Sed_level_stand_mg + Sed_exposure_days, 
##     family = binomial(link = "logit"), data = totalSS2)
## 
## Deviance Residuals: 
##    Min      1Q  Median      3Q     Max  
## -1.177   0.000   0.000   0.000   1.177  
## 
## Coefficients:
##                     Estimate Std. Error z value Pr(>|z|)
## (Intercept)         -385.217  88625.663  -0.004    0.997
## Sed_level_stand_mg     1.064    281.353   0.004    0.997
## Sed_exposure_days      4.217    968.684   0.004    0.997
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 52.3378  on 175  degrees of freedom
## Residual deviance:  5.5452  on 173  degrees of freedom
## AIC: 11.545
## 
## Number of Fisher Scoring iterations: 25

Effect estimates

## [1] 0.9886364
##    
## p     0   1
##   0 170   2
##   1   0   4
## Setting levels: control = 0, case = 1
## Setting direction: controls < cases
## Area under the curve: 0.998
##                   R2m       R2c
## theoretical 0.9996362 0.9996362
## delta       0.9966523 0.9966523

There is no significant relationship between sediment exposure and the odds of total colony mortality (concentration: GLMM z = 0.004, p = 0.997; duration: GLMM z = 0.004, p = 0.997).

Prediction figure

#Overlaying glmm results on figure, but prediction line is jagged...
pred_modSSc2 <- predict(modSSc2, newdata=totalSS2, type="response")
totalSS3 <- cbind(totalSS2, pred_modSSc2)

totalSS_plot <- ggplot(data = totalSS3) +
    geom_point(mapping = aes(
      x = Sed_level_stand_mg, 
      y = Binary_death_total,
      color = Ref)) +
    labs(x = "Sediment exposure concentration (mg/L)", 
         y = "Total mortality due to sediment exposure?",
         color = "Study") +
    scale_x_log10(limits=c(1,200),breaks=c(0.01,0.1,1,10,100,1000),
                  label=c("0.01","0.1","1","10","100","1000")) +
    geom_line(aes(x = Sed_level_stand_mg, y = pred_modSSc2), inherit.aes=FALSE) +
    theme_bw()

totalSS_plot

That plot is confusing to interpret. Let’s see if I can plot the average marginal probability, i.e., the average change in probability of the outcome across the range of the predictor of interest. This is described in some detail at the following, useful website: https://stats.idre.ucla.edu/r/dae/mixed-effects-logistic-regression/

By exposure concentration
jvalues <- with(totalSS2, seq(from = min(Sed_level_stand_mg), to = max(Sed_level_stand_mg), length.out = 100))

# calculate predicted probabilities and store in a list
pp <- lapply(jvalues, function(j) {
    totalSS2$Sed_level_stand_mg <- j
    predict(modSSc2, newdata = totalSS2, type = "response")
})

# get the means with lower and upper quartiles
plotdat <- t(sapply(pp, function(x) {
    c(M = mean(x), med = median(x), quantile(x, c(0.25, 0.75)))
}))

# add in Sed_level_stand_mg values and convert to data frame
plotdat <- as.data.frame(cbind(plotdat, jvalues))

# better names and show the first few rows
colnames(plotdat) <- c("PredictedProbability", "MedianProbability", 
                       "LowerQuantile", "UpperQuantile", "Sed_conc")
#head(plotdat)

# plot average marginal predicted probabilities
ggplot(plotdat, aes(x = Sed_conc)) + 
  geom_line(aes(y = PredictedProbability), size = 2) +
  geom_line(aes(y = MedianProbability), size = 0.5) +
  ylim(c(0, 1))

ggplot(plotdat, aes(x = Sed_conc)) + 
  geom_ribbon(aes(ymin = LowerQuantile, ymax = UpperQuantile), alpha = 0.15) +
  geom_line(aes(y = PredictedProbability), size = 2) +
  geom_line(aes(y = MedianProbability), size = 0.5) + 
  ylim(c(0, 1))

#Overlaying average marginal predicted probabilities on figure
#by Ref
totalSS_plot2 <- ggplot() +
    geom_point(data = totalSS2, mapping = aes(
      x = Sed_level_stand_mg, 
      y = Binary_death_total,
      color = Ref)) +
    labs(x = "Sediment exposure concentration (mg/L)", 
         y = "Predicted probability of total colony\nmortality due to sediment exposure",
         color = "Binary data\nby Study",
         linetype = "Predicted\nProbability") +
    scale_x_log10(limits=c(1,200), breaks=c(0.01,0.1,1,10,100,1000,loaelSS3),
                  label=c("0.01","0.1","1","10","100","1000",round(loaelSS3,digits=1))) +
    geom_ribbon(data = plotdat, aes(x = Sed_conc, y = PredictedProbability, 
                                    ymin = LowerQuantile, ymax = UpperQuantile), 
                alpha = 0.15) +
    geom_line(data = plotdat, aes(x = Sed_conc, y = PredictedProbability, 
                                  linetype = "twodash")) +
    geom_line(data = plotdat, aes(x = Sed_conc, y = MedianProbability, 
                                  linetype = "solid")) +
    geom_vline(xintercept=loaelSS3, linetype="dashed", color = "red") +
    theme_bw() +
    scale_linetype_manual(values=c("twodash", "solid"), labels = c("Median","Mean"))

totalSS_plot2

By exposure duration
jvalues <- with(totalSS2, seq(from = min(Sed_exposure_days), to = max(Sed_exposure_days), length.out = 100))

# calculate predicted probabilities and store in a list
pp <- lapply(jvalues, function(j) {
    totalSS2$Sed_exposure_days <- j
    predict(modSSc2, newdata = totalSS2, type = "response")
})

# get the means with lower and upper quartiles
plotdat <- t(sapply(pp, function(x) {
    c(M = mean(x), med = median(x), quantile(x, c(0.25, 0.75)))
}))

# add in Sed_exposure_days values and convert to data frame
plotdat <- as.data.frame(cbind(plotdat, jvalues))

# better names and show the first few rows
colnames(plotdat) <- c("PredictedProbability", "MedianProbability", 
                       "LowerQuantile", "UpperQuantile", "Sed_dur")
#head(plotdat)

# plot average marginal predicted probabilities
ggplot(plotdat, aes(x = Sed_dur)) + 
  geom_line(aes(y = PredictedProbability), size = 2) +
  geom_line(aes(y = MedianProbability), size = 0.5) +
  ylim(c(0, 1))

ggplot(plotdat, aes(x = Sed_dur)) + 
  geom_ribbon(aes(ymin = LowerQuantile, ymax = UpperQuantile), alpha = 0.15) +
  geom_line(aes(y = PredictedProbability), size = 2) +
  geom_line(aes(y = MedianProbability), size = 0.5) + 
  ylim(c(0, 1))

#Overlaying average marginal predicted probabilities on figure
#by Ref
totalSS_plot2 <- ggplot() +
    geom_point(data = totalSS2, mapping = aes(
      x = Sed_exposure_days, 
      y = Binary_death_total,
      color = Ref)) +
    labs(x = expression("Sediment exposure duration (days)"), 
         y = "Predicted probability of total colony\nmortality due to sediment exposure",
         color = "Study",
         linetype = "Predicted\nProbability") +
    scale_x_continuous(limits=c(0,max(totalSS2$Sed_exposure_days))) +
    geom_ribbon(data = plotdat, aes(x = Sed_dur, y = PredictedProbability, 
                                    ymin = LowerQuantile, ymax = UpperQuantile), 
                alpha = 0.15) +
    geom_line(data = plotdat, aes(x = Sed_dur, y = PredictedProbability, 
                                  linetype = "twodash")) +
    geom_line(data = plotdat, aes(x = Sed_dur, y = MedianProbability, 
                                  linetype = "solid")) +
    theme_bw() +
    scale_linetype_manual(values=c("twodash", "solid"), labels = c("Median","Mean"))

totalSS_plot2